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Abstract 

The evolution of Bose-Einstein condensates is amply described by the time-dependent Gross- 
Pitaevskii mean-field theory which assumes all bosons to reside in a single time-dependent one- 
particle state throughout the propagation process. In this work, we go beyond mean-field and de- 
velop an essentially-exact many-body theory for the propagation of the time-dependent Schrodinger 
equation of N interacting identical bosons. In our theory, the time-dependent many-boson wave- 
function is written as a sum of permanents assembled from orthogonal one-particle functions, or 
orbitals, where both the expansion coefficients and the permanents (orbitals) themselves are time- 
dependent and fully determined according to a standard time-dependent variational principle. By 
employing either the usual Lagrangian formulation or the Dirac-Frenkel variational principle we 
arrive at two sets of coupled equations-of-motion, one for the orbitals and one for the expansion 
coefficients. The first set comprises of first-order differential equations in time and non- linear 
integro-differential equations in position space, whereas the second set consists of first-order differ- 
ential equations with time-dependent coefficients. We call our theory multi-configurational time- 
dependent Hartree for bosons, or MCTDHB(M), where M specifies the number of time-dependent 
orbitals used to construct the permanents. Numerical implementation of the theory is reported 
and illustrative numerical examples of many-body dynamics of trapped Bose-Einstein condensates 
are provided and discussed. 
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I. INTRODUCTION 



The experimental realizations of Bose-Einstein condensates made of ultracold alkali-metal 
atoms [1-4] have stimulated a modern renaissance as to possible utilization of cold trapped 
atoms and Bose-Einstein condensates, for instance to quantum information processing [5, 
6] and interferometric-based precession measurements [7], and to the fundamental physics 
governing trapped interacting bosons [8-10]. An ultimate goal of researchers is to be able 
to design, realize, manipulate and detect desired quantum states of the many-atom system. 

Static and dynamical properties of Bose-Einstein condensates have been extensively and 
successfully explored in the community by employing the mean-field Gross-Pitaevskii theory 
[11, 12], for reviews see [8-10] and for individual applications to condensates and mixtures 
[13-19] and [20-25], respectively. Gross-Pitaevskii theory is an excellent theory for weakly- 
interacting bosons whenever a single macroscopic one-particle wavefunction is sufficient to 
describe the reality. By definition, Gross-Pitaevskii theory cannot describe phenomena such 
as fragmentation of condensates or Mott-insulator phases in optical lattices for which two, 
few or many one-particle functions are occupied. Recently, we have developed a multi-orbital 
mean-field approach to describe static and dynamical properties of fragmented condensates 
[26, 27] , thus generalizing the (one-orbital) Gross-Pitaevskii mean- field. Utilizing the multi- 
orbital mean-field has enabled us to find new phenomena associated with fragmentation, 
fermionization, quantum phases, demixing scenarios and interferences of interacting bosonic 
systems in traps and optical lattices [28-32] . 

In spite the great successes and popularity of mean-field approaches in the many-boson 
problem of degenerate quantum gases, the need and demand for practical many-body the- 
oretical approaches and computational tools is widely accepted and well documented, see, 
e.g., the review [33] and references therein. Nowadays, finite-number condensates can be 
produced (and reliably reproduced) in experiments. For instance, see the work of the Raizen 
group which demonstrated atom-number squeezing with as little as a few tens of atoms [34] , 
and the experiment of the Oberthaler group which demonstrated a single Josephson junction 
in a double- well with a few thousands of atoms [35]. Experiments as those call for delicate 
description of beyond mean-field, finite-particle-number, and many-body effects. 

The purpose of the present work is to develop an essentially-exact and numerically- 
efficient time-dependent approach for the solution of the time- dependent many-boson 
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Schrodinger equation. Recently, we have developed a general variational theory with com- 
plete self-consistency for studying stationary properties of condensates in traps on a quan- 
titative many-body level [36]. The main idea of Ref. [36] is the computation of the best 
(optimal) stationary many-body wavefunctions expressed as a linear combination of perma- 
nents. In reminiscence to our approach in the stationary many-boson problem [36], we will 
here optimize the time- dependent many-body wavefunctions according to a time-dependent 
variational principle. We term our approach multi-configurational time-dependent Hartree 
for bosons (MCTDHB). The equations-of- motion of MCTDHB have been recently posted in 
their final form in Ref. [37], where we employed MCTDHB to the popular problem of splitting 
a Bose-Einstein condensate by deforming a single well to a double-well. Applying MCTDHB 
to the ramping-up-a-barrier problem we followed the many-boson wavefunction throughout 
the splitting process and identified the role and impact of many-body excited- states on the 
splitting process. Among others, we were able to identify a new 'counter-intuitive' regime 
where the evolution of the condensate when the barrier in ramped-up sufficiently slow is not 
to the ground-state of the double-well which is a fragmented condensate, but to a low-lying 
excited-state which is a coherent condensate [37]. Here we provide the derivation of the 
MCTDHB theory, details of the numerical implementation of MCTDHB, and complemen- 
tary illustrative numerical examples. 

Evidences that employing time-dependent orbitals (permanents) in attacking the time- 
dependent many-boson Scrordinger equation beyond Gross-Pitaevskii theory is useful and 
important have already appeared in the literature. Zoller and co-workers addressed the 
ramping-up-a-barrier problem with two time- dependent orbitals of Gaussian shape whose 
positions and phase change in time [38]. More recently, Masiello and Reinhardt have used a 
time-dependent multi-configurational ansatz with two orbitals of predetermined gerade and 
ungerade symmetries to describe interferences in a symmetric double- well set-up [39]. In 
this context, we would like to stress that MCTDHB is fully variational and is not restricted 
to the number of orbitals, to a predetermined shape or symmetry of the orbitals, and to the 
geometry, dimensionality and interparticle interactions of the time-dependent many-boson 
problem. 

It is instructive to place the MCTDHB theory developed and presented in this pa- 
per in the general context of multi-configurational time-propagation approaches of many- 
particle systems. The idea to expand and optimize the time-dependent many-body wave- 
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function of distinguishable particles is long known and termed multi-configurational time- 
dependent self-consistent field approach [40-42]. A particular efficient variant of the multi- 
configurational time-dependent self-consistent field approach is the multi-configuration time- 
dependent Hartree (MCTDH) approach which has been successfully and routinely used for 
multi-dimensional dynamical systems consisting of distinguishable particles [43-46]. 

The MCTDH approach can treat efficiency and accurately static and dynamical proper- 
ties of a few-particle systems. In the latter context, we mention that very recently static 
properties of weakly- to strongly-interacting trapped few-boson systems have been studied 
on a quantitative many-body level by MCTDH [47, 48]. Yet, in treating a larger number 
of identical particles it is essential to use their statistic properties to truncate the large 
amount of redundancies of coefficients in the distinguishable-particle multi-configurational 
expansion of the MCTDH wavefunction. In this case, the challenge was first approached for 
fermionic systems where a fermionic version of MCTDH was independently developed by 
several groups [49-51], taking explicitly the antisymmetry of the many-fermion wavefunction 
to permutations of any two particles into account. 

Here we accept the respective challenge for bosons. This is in particular valuable since 
very-many bosons can reside in only a small number of orbitals. Alternatively speaking, by 
explicitly exploiting bosons' statistics it is possible to successfully and quantitatively attack 
the dynamics of a much large number of bosons with the present MCTDHB theory. A second 
importance difference in comparison with MCTDH is the nature of interparticle interactions. 
In MCTDHB we employ a general two-body interaction between identical bosons whereas 
MCTDH was designated to treat nuclear dynamics in which interactions normally involve 
several degrees-of-freedom, or coordinates. 

The structure of the paper is as follows. In section II we develop the working equations 
of MCTDHB, highlighting their appealing representation in terms of reduced one- and two- 
body densities and discussing properties of MCTDHB. In section III we present details of 
the numerical implementation of MCTDHB. Illustrative numerical examples are provided in 
section IV, whereas summary and conclusions are given in section V. Finally, complementary 
material is differed to appendices A and B. 
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II. THEORY 



The evolution of N interacting structureless bosons is governed by the time-dependent 
Schrodinger equation: 

ot N N 

H* = i-^, H( ri , r 2 ,...,r N ) = J2 K^j) + E W ( T i ~ (*) 

j=l k>j=l 

Here h — 1, rj is the coordinate of the j-th boson, h(r) = T(r) + V(r) is the one-body Hamil- 
tonian containing kinetic and potential energy terms, and W(jj — r k ) describes the pairwise 
interaction between the j-th and k-th bosons. In the most general case, the one-body poten- 
tial V(r) and the two-body interaction W(rj — r^) and, hence, the many-boson Hamiltonian 
H itself may be time-dependent. To avoid cumbersome notation in what follows, we do not 
indicate explicitly this time-dependence unless it is necessary. 

The time-dependent many-boson Schrodinger equation (1) cannot be solved exactly (an- 
alytically), except for a few specific cases only, see, e.g., Ref. [52]. Hence, theoretical and 
numerical strategies for attacking this problem are a must. 



A. Derivation of the working equations of MCTDHB 

We would like to start by constructing a formally-exact representation of the time- 
dependent many-boson wavefunction ty(t) describing the dynamics of iV identical structure- 
less bosons. To this end and following the Introduction part, we allow the bosons to occupy 
permanents made of time-dependent one-particle functions, or orbitals. Let us introduce a 
complete set of time-dependent orbitals {0fc(r, t)} which are normalized and orthogonal to 
one another at any time t, 

J 4>l{v,t)4> j {r,t)dT = 5 kj . (2) 

In what follows, it is convenient to work in second quantization formalism and introduce 
the set of bosonic annihilation operators corresponding to the orbitals {(f>k( r , t)}. This is 
conveniently done by employing the relation, 

h(t) = J ^(r,*)*(r)dr, (3) 

where *(r) = ^ k b k (t)(j) k (r,t) is the usual bosonic field operator annihilating a particle at 
position r. The bosonic annihilation and corresponding creation operators obey the usual 
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commutation relations bk{t)b-{t) — bj(t)b k (t) = 5kj at any time. Using the creation operators 
b\{t) we assemble the permanents, 

\n;t) = ^====(b\(t)) (#)) - &U) \vac), (4) 

where n = (ni,n2,n 3 , • • • ,um) represents the occupations of the orbitals that preserve the 
total number of particles n\ + n 2 + n 3 + ■ ■ ■ + n M = N, M is a number of the one-particle 
functions, and \vac) is the vacuum. 

In the multi-configuration time-dependent Hartree approach for bosons (MCTDHB) the 
ansatz for the many-body wavefunction \&(t) is taken as a linear combination of time- 
dependent permanents (4) [37], 

\y(t)) = J2Cn(t)\n;t), (5) 

n 

where the summation runs over all possible configurations whose occupations n preserve the 
total number of bosons N. Of course, in the limit M — > oo the set of permanents {\n;t)} 
spans the complete iV-boson Hilbert space and thus expansion (5) is exact. So where is 
the advantage of utilizing an expansion with time- dependent permanents? In practice, we 
have of course to limit the size of the Hilbert space exploited in computations. Now, by 
allowing also the permanents to be time-dependent we can use much shorter expansions 
than if only the expansion coefficients are taken to be time- dependent, thus leading to a 
significant computational advantage. We would like to highlight that in representation (5) 
both the expansion coefficients {Cft(t)} and orbitals (0fc(r,t)} comprising the permanents 
|n; t) are independent parameters. To solve for the time-dependent wavefunction ^(t) means 
to determine the evolution of the coefficients {C^(t)} and orbitals {0fe(r, t)} in time. 

To derive the equations-of motion governing the evolution of {C^(t)} and {(f>k(r,t)} we 
need to employ a time-dependent variational principle. Two such variational principles are 
utilized here, both leading to the same result, of course, but highlighting complementary 
aspects when treating the time-dependent many-boson problem. Specifically, in the bulk of 
the paper below we employ the usual Lagrangian formulation [53, 54], whereas how to derive 
the MCTDHB equations-of-motion starting from the Dirac-Frenkel variational principle [55, 
56] is deferred to appendix B. The main difference between the two (equivalent) formulations 
employed here is that, in the Lagrangian formulation we are to take first expectation values 
and only subsequently perform the variation which somewhat simplifies the algebra, whereas 
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in the Dirac-Frenkel formulation the situation reverses, i.e., variation of ty(t) precedes the 
computation of matrix elements. 

In the framework of the Lagrangian formulation [53, 54], we substitute the many-body 
ansatz (5) for into the functional action of the time-dependent Schrodinger equation 
which reads: 

(6) 



S[{Cn(t)},{<f> k (*,t)}} 




k,j 

The time-dependent Lagrange multipliers fikjit) are introduced to ensure that the time- 
dependent orbitals 4>k( r , t) remain orthonormal throughout the propagation, see (2) and 
appendix B. The next step is to require stationarity of the action with respect to its 
arguments {Ca(t)} and {(j) k (r,t)}. This variation is performed below separately for the 
coefficients and for the orbitals, recalling that they are independent variational parameters. 



1. Variation with respect to the expansion coefficients {C^(i)} 



To take the variation of the functional action (6) with respect to the expansion coefficients 
{Cn(t)} we first express the expectation value in the action in a form which explicitly depends 
on the expansion coefficients, 

.dC, 



n'; t ) Cft> — i- 



dt 



(7) 



It is now straightforward to perform this variation which gives, 



dS[{C H (t)},{Mr,t)}} 







E 



n; t 



H-i 



d_ 

dt 



dcm 

Defining the time-dependent matrix 7i(t) the elements of which are 

d 



n ; t ) Cfi 



. dC H 
' dt ■ 



(8) 



7~Lnri' (t) 



n; t 



H-i 



dt 



n ; t 



(9) 



and collecting the expansion coefficients Cft{t) in the column vector C(t), Eq. (8) can be 
written in a compact matrix form as: 

.dC(t) 



H(t)C(t) 



dt 



(10) 



Let us discuss the properties of Eq. (10). The equation-of-motion (10) is a first-order differ- 
ential equation in time. The matrix 7i(t) multiplying the vector C(t) on the left-hand-side 



is time-dependent, since it is evaluated with time-dependent permanents \n;t) and \n';t), 
comprised themselves of time-dependent orbitals. 

In the course of evolution, the many-body wavefunction \I/(t) should, of course, re- 
main normalized since for self-adjoined Hamiltonians H the evolution is unitary. \^{t)) = 
Xln Cft(t) \n; t) would remain normalized if the vector of coefficients C(t) remains normalized 
in time throughout propagation via Eq. (10) (the orbitals comprising the permanents remain 
orthonormal to one another by virtue of the Lagrange multipliers, see subsection II A 2 for 
more details). For this, the matrix 7i(t) should be hermitian. 

Below, we first evaluate the matrix 'Hit) explicitly and then show that it is always her- 
mitian. Derivation of the rules to evaluate matrix elements of one- and two-body operators 
between two general permanents can be found in [36]. Here the matrix elements with re- 
spect to two general time- dependent permanents \n; t) and \n';t) are displayed in their final 
form. Noticing that in second quantization the time-derivative can be treated as a one-body 
operator, 

(11) 



the non-vanishing matrix elements of the one-body operator h 



■th 



are given by 



ft: t 



n k -t 



h - i at 



n: t 



ft; t 



M 

i=i 



mi — ^tt 
dt 



sj (n k + l)n q 



h 



kq 



dt) 



kq 



, k ^ q, 



(12) 



where the permanent denoted by \n^\ t) differs from \n; t) by an excitation of one boson from 
the g-th to the fc-th orbital. The corresponding matrix elements of the two-body operator 
W between two permanents \n; t) and \n'\ t) are collected for convenience in appendix A. 
The time-dependent matrix elements of the one- and two-body operators with respect to 
orbitals needed in this work are given by: 

h kq = J <j>l(T,t)h(r)<i> q (T,t)dT, 

W ksql = j j (j>l(r,t)(j>* a (T f ,t)W(T -I f )(j> q (T,t)(l> l (l f ,t)dTdT f , 
W ks [ql] = W ksq l + W ks l q . (13) 



Note the plus sign in W ks [ q i] which is due to bosons' statistics. 

With the elements of lH.it) explicitly expressed, we can return now to the question of its 
hermiticity. Clearly, the contribution originating from the matrix elements of the Hamil- 
tonian, hk q and Wk sq i, is hermitian for any set of of functions {0&(r, t)} belonging to the 
domain of the Hamiltonian. To show that the matrix representation of is also hermitian 
we make use of the fact that the orbitals {</>fc(r, t)} are kept normalized and orthogonal to 
one another at any time. Taking the time derivative of Eq. (2) we readily get 

. d 



j (f>l(T,t)<f> q (T,t)dT 



Summarizing, we have shown that the matrix Hit) is hermitian for a general set of orbitals 
{0fc( r , t)} that are kept normalized and orthogonal to one another at any time. Consequently, 
the evolution of the expansion coefficients {Cft(t)} is always unitary, namely, that an initially 
normalized expansion coefficient vector C(t) propagated by Eq. (10) remains normalized at 
any time. This result, together with Eq. (2), guarantees that an initially normalized many- 
body state ^(O) remains normalized at any time. 



2. Variation with respect to the orbitals {4>k{ Y -> t)} 



In the present subsection we derive the equations-of-motion governing the evolution of the 
orbitals {4> k (r,t)}. Now, it is helpful to express the expectation value of H — appearing 
in the functional action (6) in a form which allows one for a direct functional differentiation 
with respect to 0jt(r, t). 

To this end, we write the operator H — in second quantization form: 



H-i 



d_ 

dt 



k,q 



hkq — 







dt) 



kq 



+ 



\ E b lb%biW ksql , 



(15) 



k,s,q,l 



where the matrix elements of the one- and two-body operators are given in Eq. (13). In 
calculating the expectation value of (15) with respect to the many-boson wavefunction ty(t) 
it is gratifying to make use of the one- and two-body density matrices Pk q (t) — \ \P 



h'h 



and pksqi{t) = (* 



blb\b q k 



respectively. Given the (normalized) wavefunction ^{t), the 



9 



one-body density matrix is given by 



p(ri|ri;£) = N J ^*(r[, r 2 , . . . , r^; t)tf(ri, r 2 , . . . , rjv; t)dr 2 dr 3 ■ ■ ■ dr N 



M 



fc,g=l 



where the matrix elements of the density explicitly read 



Pkk(t) = y]CgCjtWfc, 

n 



(16) 



(17) 



It is convenient to collect these matrix elements as p(t) = Similarly, the two-body 

density matrix associated with ty(t) is given by 

p(r u r 2 |ri, r 2 ; f) = iV(iV - 1) ^ r' 2 , r 3 , . . . , rjv ; f)tf (n, r 2 , r 3 , . . . , r N ; t)dr 3 ■■■dr N = 



M 



fe,s,g,i=l 

(18) 



where the matrix elements pk sq i are collected for convenience in appendix A. 

Combining the above ingredients we obtain for the expectation value of the operator 
H — ijr appearing in the factional action (6), 



M 
k,q=l 



^kq 



( d \ 1 1 M 

^ ' fc <?J k,s,q,l=l fx 



•w- (19) 



Note the appealing appearance of representation (19) in which the only explicit dependence 
on the orbitals {0fc(r, t)} is grouped into the matrix elements hk q , (^§i) kq an d Wk sq i of the 
one- and two-body operators h — i^- t and W , respectively, whereas the elements pk q and pk sq i 
of the density matrices do not depend explicitly on the orbitals. Consequently, it is now 
straightforward to perform variation of the functional action (6) with respect to the orbitals 
{0fc(r, t)} which gives the following set of coupled equations-of- motion, k — 1, . . . , M, 



SS[{C H (t)},{Mr,t)}] 
6<f>l(v,t) 



= 



M 

E 

q=l 



d 



M 



Pk q [ h - i—J + X Pks q iW s i 

' s,l=l 



M 



10?) = J2pkj\<i>j), 

(20) 
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where 



(21) 



are local time-dependent potentials. 

To proceed, we can use the constraints (2) that the <f>k(r,t) are functions ortho normal to 
one another in order to eliminate the Lagrange multipliers /%■(£) from Eq. (20). By taking 
the scalar product of {<fij\ with (20), the resulting /%(t) take on the form, 



M 



9=1 



h j9 ~ ( ^ 



:i<i 



M 



+ PksqlWjsgl 
s,l=l 



(22) 



Substituting Eq. (22) into (20), employing the identities 



M 

E 

9=1 



/ Q\ M 

Pkq ( ^ - i qT J + ^ PksqlW sl 
^ ' s,l=l 



M 



M \ M 

i'=l / 9=1 



10?) |$) 
J=l 
M 



Pkq\h-i 



s,/=l 



!</.,), fc = l,...,M, (23) 



and multiplying from the left with the inverse of the one-body density and summing over 
^2 k {p(t)}~k, we arrive immediately at the following form of the equations-of-motion of the 
orbitals (f>j(r, t), j = 1, . . . , M: 

M 

£ {p{t))]k pks q iW s i\4> q ) 



Pi 



fa 



k,s,q,l=l 



M 



i'=i 



(24) 



Here and hereafter we use the shorthand notation 3 - = -J^.. Examining Eq. (24) we see 
that eliminating the Lagrange multipliers Hkjif) has emerged as a projection operator P 
onto the subspace orthogonal to that spanned by the {0fc(r, £)}. This projection appears 
both on the left- and right-hand-sides of (24), making (24) a cumbersome coupled system 
of integro-differential non-linear equations. 

Fortunately, due to the invariance properties of the ansatz (5), see discussion in subsection 
II B below, we can always perform a unitary transformation without introducing further 
constraints on the orbitals such that [43, 44] 



k\<t>q) = 0, k,q=l, 



,M, 



(25) 
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are satisfied at any time. Obviously, if conditions (25) are satisfied at any time, the or- 
thogonality constraints (2) are also satisfied. This representation simplifies considerably the 
equations-of-motion (24) for the orbitals (f)j(r, t), j — 1, . . . , M: 

M 

h\<f>,)+ Yl {p^yjk pk Sq iw sl \4> q ) 



k,s,q,l=l 



M 



3>=i 



(26) 



The P remaining on the right-hand-side of Eq. (26) makes it clear that conditions (25) 
are indeed satisfied at any time throughout the propagation of {4>k(r,t)}. In practice, 
the meaning of these conditions is that the temporal changes of the {0 fe (r, t)} are always 
orthogonal to the {0fe(r,t)} themselves. This property also utilized in MCTDH [43, 44] 
generally makes the time propagation of Eq. (26) robust and stable and can thus be exploited 
to maintain accurate propagation results at lower computational costs. Additionally with 
conditions (25), Eq. (10) now reads: 



H(t)C(t) 



.dC(t) 



H 



(27) 



n'-t). 



where H(£) is the Hamiltonian matrix the elements of which are Hrm'(t) = \n;t 
The coupled equation sets (26) for the orbitals {(/>,,■ (r, £)} and (27) for the expansion coeffi- 
cients {Cfi(t)} are at the heart of the multi-configurational time-dependent Hartree theory 
for bosons - a formally-exact and practical representation of the time-dependent many-boson 
Schrodinger equation (1). As mentioned above, we term our theory in short MCTDHB(M), 
where M stands for the number of time-dependent orbitals comprising the permanents. Of 
course, in the limit M — > oo the set of permanents (|n;t)} spans the complete iV-boson 
Hilbert space and thus ty(t) is exact. In fact, inspecting Eq. (26) and the projector P 
therein tells us that in this limit the orbitals {</>fc(r,t)} become time- independent, similarly 
to the situation in the general MCTDH theory [43, 44]. So where is the advantage of uti- 
lizing time-dependent permanents? In practice, we have of course to limit the size of the 
Hilbert space used in computations. By allowing the permanents to be time-dependent, 
variationally-optimized quantities we can use much shorter expansions than if permanents 
comprised of fixed-shape orbitals are employed, thus leading to a significant computational 
advantage. 
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B. Properties of MCTDHB and its working equations 



As mentioned in the Introduction, it is gratifying to note that the present many-body 
propagation theory - MCTDHB - adapts to identical bosons the multi-configurational time- 
dependent Hartree (MCTDH) approach routinely used for multi-dimensional dynamical sys- 
tems consisting of distinguishable particles [43-46] . By explicitly exploiting bosons' statistics 
and the fact that only a two-body interaction is needed it is possible to successfully and quan- 
titatively attack the dynamics of a large number of bosons with MCTDHB. Hence, many of 
the properties of MCTDHB are inherited from MCTDH. In this section we concentrate on 
and expand those properties that are required for our needs. 

The ansatz (5) for the many-boson wavefunction includes all possible permanents \n;t) 
assembled when distributing N bosons over M (time-dependent) orbitals <f>k(t). Since the 
above defines complete Hilbert subspaces, \^{t)) possesses invariance properties, i.e., it does 
not have a unique representation in terms of the orbitals and coefficients. Specifically, we 
can introduce an M x M unitary matrix U(£) = {Uk q (t)}, define a new set of orthonormal 

n:t) with them. 



orbitals (j> q (t) = Uk q (t)(/)k(t) : and assemble all possible permanents 
Correspondingly, the vector of coefficients is transformed and (5) can be rewritten to express 
this invariance, 



The size of the Hilbert subspace remains the same, of course. To express that we use the 
same vector of enumeration n. 

The invariance of \^{t)) to unitary transformations of the orbitals has been used in sub- 
section II A 2 to simplify the form (24) and (10) of the equations-of-motion to the respective 
form (26) and (27) of MCTDHB which is amenable to numerical implementation. This has 
been achieved by employing the differential conditions (25). We stress that an invariance 
like (28) has been used by the developers of MCTDH in order to choose the differential 
form (25) to simplify the MCTDH equations-of-motion [43, 44]. Although known from their 
original MCTDH papers for distinguishable particles [43, 44], it is deductive to stipulate 
an explicit proof that employing condition (25) throughout the time evolution does not 
introduce further constraints or approximations on the variational treatment of the time- 
dependent many-boson (many-body) Schrodinger equation. To this end, we show by direct 
construction that there is a unitary matrix U(t) that carries a general set of time-dependent 
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orthogonal orbitals to a specific set of time-dependent orthogonal orbitals satisfying Eq. (25). 

To obtain an equation-of- motion for U(t), we start from the hermitian matrix D kq (t) = 
i (^(fik(t)\<fiq(t)^ , whose elements are not necessarily zero numbers, and wish to com- 
pute the unitary transformation <f> q (t) = J^fcli U kq (t)(j) k {t) which guarantees the rela- 
tions D kq (t) = i (^(f) k (t)\(f) q (t)^ = 0. For this, we substitute the time-derivative 4> q {t) = 

Y^k=i \Uk q {t)<t>k{t) + U kq (t)<j) k (t)^ into the equation D kq {t) = and solve for the matrix 
U(£) which ensures these conditions at all times. Making explicit use of the assumption 
that U(t) is unitary we immediately get and symbolically integrate, 

M 

iU sq (t) ^-J2 D ^(t)Ukq(t) =► U(f) = e + ^ D ^'U(0). (29) 
fc=i 

Obviously, since D(t) is an hermitian matrix, U(t) remains unitary at all times if and only 
if the initial condition U(0) is an unitary matrix. We show now by direct construction that 
U(0) = lim T ^ U(r) is unitary and unique. For this we diagonalize the hermitian matrix 
D(£) with the help of the unitary matrix T(£), 

Tt(f)D(f)T(f)=d(f), (30) 

where d(t) is the diagonal matrix of the eigenvalues of D(t). From this it is not hard to see 
that in the limes r — > 0, 

M 

U(r) = T(0)e + - d (°), J ff (r) = ^ U kq {r)4> k {r) =► D kq (r) = 0, (31) 

k=l 

which concludes our proof. 

Thus, for any set of initial conditions {4>k(0)} we can choose the constraint (25), i.e., 
D kq (t) = 0, and safely work with equations-of- motion (26) and (27) of MCTDHB without 
loss of generality. In fact, as has been shown in the original MCTDH work [44], we can use 
a larger class of constraints, namely D kq (t) = (<f) k \g\ 4> q ), where g is a self-adjoint operator 
in the subspace of orbitals {(f) k (r,t)}. The proof that such a choice only amounts to a time- 
dependent unitary transformation of the orbitals and leads to no further restrictions follows 
exactly the same lines as above. Here we just write the final result for the equations-of- 
motion which take the form, j — 1, . . . , M: 

M 

^)=g\ ( p j ) + P [ ] 

k,s,q,l=l 

W(t)C(f)=i^W (32) 
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where now TL{t) is the time-dependent matrix with elements TCnn'(t) — (n;t 



H-g 



n':t). 



The next property we would like to address is that of the MCTDHB energy, 
ty(t) H ty(t)^. For time-independent Hamiltonians, as has been shown for MCTDH it- 
self [43, 44], the energy is constant in time as it should be for time-evolution with such 
Hamiltonians. For a time-dependent Hamiltonian H(t), it is not difficult to show by direct 
differentiation of the MCTDHB energy and utilizing without loss of generality relations (25) 
that, 



H 



dH 



dt 



(33) 



This appealing relation can be used in numerical calculations to monitor the degree of 
accuracy of integration. 

Finally, there are two points worth mentioning. First, we can also propagate the MCT- 
DHB equations in imaginary time and compute for static (time-independent) traps self- 
consistent ground and excited eigenstates of bosonic systems, since in that case MCTDHB 
boils down to the general variational many-body theory for interacting bosons (MCHB) 
developed recently in [36]; see [57] for the corresponding MCTDH case. Second, with an 
essentially-exact and full many-body theory for the dynamics of bosonic systems we can now 
investigate strategies for approximate many-body time-dependent theories, e.g., when not 
all coefficients {Cft(t)} are employed in the ansatz for ^(t). 



III. NUMERICAL IMPLEMENTATION 



In the present section we describe the numerical implementation the developed MCTDHB 
theory for systems of cold bosonic atoms. Here we would like to mention that in our 
independent technical implementation of MCTDHB we have borrowed much of the ideology 
of the implementation of the general MCTDH theory [43-46]. Still, we reiterate two major 
differences, being the direct utilization of bosons' statistics and exploitation of two-body 
interactions which translate to the appearance of the two-body density in the MCTDHB 
theory. 

To integrate the time-dependent many-boson Schrodinger equation means to find the 
many-body wavefunction at time t by specifying the initial conditions, i.e., the many-body 
wavefunction at time zero. In the framework of the MCTDHB theory we have to solve 
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simultaneously Eqs. (26) and (27) by specifying the initial set of expansion coefficients C(0) 
and respective orbitals <fr(0) at t — 0. Here and hereafter we use the column vector <f> = {<f) k } 
to group together the orbitals. In reality the choice of the initial guess is quite a complicated 
task which depends on a specific experimental setup and/or on the experimental sequence 
applied to the initial atomic cloud. Typically as an initial guess one uses a many-body 
solution of the stationary many-boson Schrodinger equation. 

According to the MCTDHB derivation, Eq. (27) determines the evolution of the expan- 
sion coefficients for a given set of the orbitals and Eq. (26) governs the evolution of the 
orbitals for a given set of the expansion coefficients. However, it is important to note that 
Eq. (27) determining the evolution of the expansion coefficients does not depend explicitly 
on the orbitals, rather it depends on the one- and two-body matrix elements h kq and Wk sq i, 
see (13). Analogously, the expansion coefficients enter Eq. (26) only implicitly via the el- 
ements of the one- and two-body densities p kq and Pksqi, see (16) and (18). Let us display 
the MCTDHB equations (26) and (27) in a form where these functional dependencies are 
explicitly indicated: 

ic(t) = o 1 {h kq m)] , w ksql [<t>{t)\ , c(*)} 

i<f>(t) = 2 { Pkq [C(t)\ , Pksql [C(t)\ , <f>(t)} , (34) 

where the quantities Oi and 2 represent column vectors with functional dependences as 
indicated. Before prescribing how to solve this coupled system simultaneously, let us first 
analyze each of these equations. The first of these equations is linear (see Eq. (27)) with 
respect to expansion coefficients C and thus, if the matrix elements h kq and W ksq i are 
given, can be effectively solved by integrators explicitly designed for linear equations, such 
as short iterative Lanczos (SIL) integrator [58]. The second equation for the orbitals is 
more sophisticated, because apart of the differential T and local V and W st (see Eq. (21)) 
operators, it also contains the projection P classifying it as an integro-differential equation. 
Again, if the elements of the one- and two-body densities p kq and p ksq i are given, this non- 
linear integro-differential equation can be propagated by means of general variable-order 
integrators, such as Adams-Bashforth-Moulton (ABM) predictor-corrector integrator [59] 
which is our choice here. Now we are ready to integrate simultaneously the coupled system 
(34). We recall that these equations are coupled through the time-dependent quantities 
h kq {t) and W ksq i(t) which depend only on the orbitals, and the quantities p kq (t) and p ksq i(t) 
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which depend only on the expansion coefficients. Therefore, if t is discretized these numbers 
can be kept unchanged during each time-step. Going along the MCTDH lines we apply 
a second order integration scheme with estimation on discretization error and adjustable 
time-step size r [45, 46]. To propagate the coupled system (34) from to r we have used 
the following flowchart: 

• Having at hand the initial conditions, i.e., the set of orbitals 0(0) and expansion 
coefficients C(0) at t — 0, one computes the quantities h kq (0), Wk sq i(0) and p kq (0), 
Pksgi(0), respectively. 

• Propagate C(0) -> C(r/2) with SIL using h kq {0) and W kaql (0). 
Evaluate p kq {r/2) and Pk sq i(r/2) using C(r/2). 

• Propagate 0(0) — > 4>{t/2) with ABM using p kq (r/2) and p ksq i(T/2). 

• For error estimation; Propagate 0(0) — > 0(r/2) with ABM using p kq (0) and p ksq i(0). 

• Propagate 0(r/2) -> 0(r) with ABM using p kq (r/2), p ksql {r/2). 
Evaluate h kq (r) and W ksq i{r) using 0(r). 

• Propagate C(r/2) — > C(r) with SIL using h kq {r) and H^ ag i(r). 

• For error estimation; Backward propagate C(0) <— C(r/2) using h kq (r) and H^ sg j(r). 

The differences C(0) — C(0) and 0(r/2) — 0(r/2) are used to estimate the discretization 
error and to adjust the time-step size r. 

We have so far fully implemented MCTDHB(M = 2), i.e., with two orbitals, in one- 
dimension, and for a general two-body interparticle interaction W(x — x') which is con- 
veniently represented in a product form, W[x — x') = J2i u i{ x ) u i( x ')- ^ n ^ ne context 
of quantum gases we have also implemented the popular contact interparticle interac- 
tion W(x — x') = \q5(x — x'), see Refs. [8-10]. To represent the one-body Hamiltonian 
h(x) = T(x) + V(x) and the orbitals we employ the discrete variable representation (DVR) 
technique [60] based on harmonic oscillator, sinusoidal or exponential functions. In the DVR 
basis the kinetic-energy differential operator T{x) corresponds to a non-diagonal matrix, the 
potential V(x) is diagonal, and the orbitals are represented by column vectors on the re- 
spective DVR grid. Therefore, in the DVR technique a differentiation of the wavefunction 
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is equivalent to matrix to vector multiplication, whereas its integration corresponds to sum- 
mation of the elements of the respective vector(s) with the corresponding DVR weights [60]. 
Finally, the accuracy of computations and respective integration efficiency of the numerical 
schemes described above depend on the specific many-boson problem under investigation, 
see for specific examples the following section IV. 

IV. ILLUSTRATIVE NUMERICAL EXAMPLES 

The first application of MCTDHB is for the dynamics of splitting a condensate when 
ramping-up a barrier such that a double-well is formed and can be found elsewhere [37]. 
Before turning to the examples of this work we first briefly summarize the main outcome 
of our application in [37]. We have found that the dynamics of splitting when ramping-up 
a barrier depends on the duration of the process and on the (effective) interaction strength 
between the bosons. There are (at least) two distinct regimes: (i) an adiabatic regime 
where the initial condensed ground-state evolves towards the ground two-fold fragmented 
eigenstate of the final double- well potential and asymptotically approaches it with increasing 
rumping-up time; (ii) an inverse regime where the initial condensed state evolves towards 
the ground two-fold fragmented eigenstate only for short rumping-up times, while for slow 
rumping-up processes the time-dependent state stays condensed during all the evolution and 
thereby evolves to a non-ground many-body eigenstate, see [37] for details. 

The main purpose of our studies performed here is to compare the mean-field dynamics 
calculated by the time-dependent Gross-Pitaevskii equation and the many-body dynamics 
calculated with MCTDHB, and thereby study some differences between the corresponding 
dynamics in simple and deductive examples. 

In the numerical examples below we consider a many-boson system in one dimension. We 
choose for convenience a length scale L such that the energy unit is = 1, where m is the 
mass of a boson, and then translate the time-dependent many-boson Schrodinger equation 
to dimensionless units. The one-body Hamiltonian then reads h(x) = + V(x). Of 

course we have also to choose a specific shape for the interparticle interaction and we do so 
by taking the popular contact interaction W[x — x') = \ 5(x — x f ), see Refs. [8-10]. The 
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resulting two-body matrix elements (13) and time-dependent potentials (21) simplify, 



W ksq i(t) = A J (j)* k (x,t)(f)* s (x,t)(f) q (x,t)(f)i(x,t)dx, W ks[q i] = 2W ksqh 
W sl (x,t)=\ (i>* s (x,t) ( f> l (x,t). (35) 



With these quantities computed, Eq. (26) for the evolution of the orbitals reads, 

i 4) = p 



M 

j'=i 



M 

h \<t>j) + A {til* Pksql^l \<f>q) 

k,s,q,l=l 

and the matrix elements of the Hamiltonian if^/(i) = (n\t H n';tj needed for the prop- 
agation of the coefficients in (27) are readily evaluated. 

As illustrative numerical examples we consider the following scenario. We prepare 
a system comprising of iV bosons in the ground-state \&(0) of the double- well potential 
Vq(x) = + 8e~ x2 ^ 2a2 \ a = 2.6. The initial state \&(0) is computed by imaginary time 
propagation of MCTDHB(2) and separately by the Gross-Pitaevskii equation for the sake 
of later comparison. At time t = we halve the barrier between the two wells and ad- 
ditionally translate the whole potential to the left. The resulting double- well potential is 
V(x) = ^Hl! _|_ 4 e -(z+2) 2 /(2<x 2 ) ? or = 2.6. The many-boson state ^(0) in which the system 
is prepared is obviously not in a stationary state any more and the interacting system is let 
to evolve in time. The time-dependent many-boson wavefunction is respectively computed 
by now real time propagating of MCTDHB(2) and of the Gross-Pitaevskii equation. Two 
systems are considered. The first with N = 100 bosons and Ao = 9.99/99 = 0.01009, the 
second with iV = 1000 bosons and A = 9.99/999 = 0.0100. Note that the two systems are 
characterized by the same 'mean- field' factor \ (N — 1) = 9.99. This means that both sys- 
tems would show identical mean-field dynamics, since the only factor concerning the number 
of bosons and their mutual interaction entering the Gross-Pitaevskii theory is the product 
Xo(N-l) [8-10]. 

Before moving to present and analyze the results, we would like to record some technical 
data used in our calculations. For the DVR we use 257 sinusoidal functions, which we 
find to be sufficient to fully converge the present results. The ABM integrator employed to 
propagate the orbitals (see section III) is set to 7-th order. For N = 100 bosons the maximal 
size of the SIL subspace needed for convergence is as low as 8 vectors, whereas for iV = 1000 
bosons it comprises of only 25 vectors! Average values for the integration time-step size r are 
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0.007 and 0.002 for iV = 100 and N = 1000 bosons, respectively. Finally, in all calculations 
and with the above parameters the integration error of ty(t) is 10~ 8 , combining the errors 
introduced by the SIL and ABM integrators (see section III). 

We compare the time-dependent dynamics of the mean-field Gross-Pitaevskii and of 
the many-body MCTDHB(2) approach. Since visualization of the time-dependent many- 
body wavefunction ty(t) is quite cumbersome, we plot in Fig. 1 snapshot of the density, 
p(x,t) = q =\ Pfcg(^)0fc(^, at different times and compare with the respective 

Gross-Pitaevskii theory. Furthermore, in Fig. 2 we display the natural occupation num- 
bers Pj(t), j = 1,2 of i.e., the eigenvalues of the corresponding reduced one-particle 
density, p(x\x';t) = X^=i Pj^'t'^ ( x ' ■> t)<j>f°(x, t), at each point in time. Of course, in 
Gross-Pitaevskii theory there is only one occupation number, p\ = 100%. 

Let us analyze the densities shown in Fig. 1. Since the barrier between the two potential 
wells is initially quite high, the corresponding three densities coincide at t — 0, see top panel 
of Fig. 1. At t — 3.0 we already see an observable difference. The Gross-Pitaevskii density 
acquires wiggles due to collisions with the potential walls. On the other hand, the densities 
in the many-body calculations show almost no wiggles for N = 100 and N = 1000 bosons. 
By t = 50.0 all three densities have substantially distorted from the initial conditions and 
exhibit multiple density wiggles of various sizes and depths. All densities are different from 
one another. In particular, the systems with N = 100 and N = 1000 bosons which, as we 
have mentioned above, have the same time evolution on the mean-field level are actually 
showing quite distinct evolutions on the many-body level, see Fig. 1. 

Next, let us analyze the natural occupation numbers of ^(t). In the Gross-Pitaevskii 
dynamics the initial condition, which is the ground state of the potential Vq(x) computed 
by imaginary-time propagation of the Gross-Pitaevskii equation, is fully coherent and hence 
there is only one occupation number, p\ = 100%. This property of the condensate obviously 
does not change in time when the time-dependent Gross-Pitaevskii equation governs the 
dynamics. Contrastly, on the many-body level the situation differs substantially. First, the 
initial conditions for N = 100 and N = 1000 bosons, which are the respective ground states 
of the potential Vq(x) computed by imaginary-time propagation of the MCTDHB(2) equa- 
tions, actually correspond to fragmented condensates, see Fig. 2 at t — 0. Most importantly, 
in the course of the many-body time evolution of the bosonic systems the natural-orbital 
occupations change in time by at least 20%. This suggests an intricate and important prop- 
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erty of the many-body dynamics which by definition cannot be accessed by the mean-field 
dynamics. 

All in one, the numerical results of the application of MCTDHB(2) to simple trapped 
bosonic systems demonstrate that the many-boson dynamics in traps is manageable and 
intriguing, and that it can also be quite different from the mean-field dynamics. 

V. SUMMARY AND CONCLUSIONS 

The evolution of Bose-Einstein condensates has been extensively studied by the well- 
known time-dependent Gross-Pitaevskii equation which assumes all bosons to reside in the 
same time-dependent orbital and hence the condensate to be coherent at all times. In 
this work we address the evolution of condensates on the many-body level, and develop and 
report on an essentially-exact and numerically-efficient approach for the solution of the time- 
dependent many-boson Schrodinger equation. We term our approach multi-configurational 
time-dependent Hartree for bosons (MCTDHB). 

In the MCTDHB, the ansatz for the many-boson wavefunction is taken as a linear com- 
bination of all possible time- dependent permanents made by distributing N bosons over 
M orthogonal time- dependent orbitals. The evolution of the many-body wavefunction is 
then determined by utilizing a standard time-dependent variational principle - either the 
Lagrangian formulation or the Dirac-Frenkel variational principle. Performing the varia- 
tion, we arrive at two sets of coupled equations-of-motion, one for the orbitals and one for 
the expansion coefficients. The first set comprises of first-order differential equations in 
time and non-linear integro-differential equations in position space. The second set consists 
of first-order differential equations with time-dependent coefficients. MCTDHB naturally 
relates - by performing propagation in imaginary time - to the recently developed and suc- 
cessfully employed general variational mean-body theory with complete self-consistency for 
stationary many-boson systems (MCHB). From another instructive perspective, MCTDHB 
can be seen as specification of the successful multi-dimensional wavepacket-propagation ap- 
proach MCTDH to systems of identical bosons. By explicitly exploiting bosons' statistics 
and the fact that only a two-body interaction is needed it is now possible to successfully and 
quantitatively attack the dynamics of a much larger number of bosons with the developed 
MCTDHB theory. 
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With an essentially-exact and full many-body theory for the dynamics of bosonic systems, 
we can now investigate the many-body correlation effects coming atop the time-dependent 
Gross- Pit aevskii [11, 12] and multi-orbital [27] mean-field dynamics of unfragmented and 
fragmented Bose-Einstein condensates, respectively. Another promising direction would be 
to utilize MCTDHB for devising and investigating strategies for approximate many-body 
time-dependent approaches such as Hilbert-space truncation schemes. 

Finally, illustrative numerical examples of the many-body evolution of condensates in 
a one-dimensional double-well trap are provided, demonstrating an intricate dynamics of 
the density and natural orbitals as time progresses which are not accounted for by the 
respective Gross-Pitaevskii mean-field dynamics. This is the tip of the iceberg of many-body 
dynamical properties of Bose-Einstein condensates and interacting bosons which await many 
more investigations with MCTDHB. 
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APPENDIX A: MATRIX ELEMENTS OF TWO-BODY OPERATORS AND OF 
THE TWO-BODY DENSITY MATRIX 

Matrix elements of the two-body operator W between two general permanents \n; t) and 
\n';t). The permanent denoted by \n^;t) differs from \n;t) by an excitation of one boson 
from the q-th to the k-ih orbital, whereas |^f;t) differs from \n;t) by excitations of two 
bosons, one boson from the q-th to the fc-th orbital and a second boson from the Z-th to the 
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s-th orbital: 
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n k + l)(n k + 2)(n q - l)n q W kkqq 



\J {n k + l)(n k + 2)n q ni W kkq i 

yj {n k + l)(n s + l)(n q - l)n q W ksqq 

sj {n k + \){n s + l)n q m W ks M . 



(Al) 



Different indices k, s, q, I cannot have the same value. All other matrix elements vanish. 
Matrix elements of the reduced two-body density p(r 1 ,r 2 \r' 1 ,r' 2 ;t) 
V(t) *Vi)*V 2 )*( r i)*( r 2)| of the many-boson state |tf(t)) = £ B C s (f) \n;t): 



Pkkkk — 53 C* H Cfi (n\ — n k ) , 

n 

Pksks — 53 CfjCfjn k n s , 



Pkkqq = 53 C n C ^kl\l ^ ~ l ) n k( n i + !)K + 2 )> 
n 

Pfcfcfc/ = ^QQ t K - l)\Jn k {ni + 1), 

n 

= 53 C r-t C n%n s ^/ n k (n s + 1), 



PfcM = 53 V (wfc - l)n-k(n q + l)(nj + 1), 



Pksqq = 53 ^Cn^ \J n k n s {jl q + l)(n q + 2), 
n 

Pfcss/ = 53 C H C a l k n a y/n k (ni + 1), 

n 

Pfe S g« = 53 C « C «2 , (1 \f^nJn g + + 1). 



(A2) 



Different indices fc, s, g, Z cannot have the same value. All other matrix elements can be 
computed due to the symmetries of the two-body operator, p ksq i = p ks i q = p sk i q = p skq i, and 
its hermiticity, p* ksql = p qlks . 
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APPENDIX B: DERIVATION OF THE MCTDHB EQUATIONS-OF-MOTION 
STARTING FROM THE DIRAC-FRENKEL VARIATIONAL PRINCIPLE 



The Dirac-Frenekl variational principle reads [55, 56]: 



6*(t) 



d 



*(0) = o. 



(Bl) 



Given the MCTDHB ansatz for the many-boson wavefunction, = Y^nCn(t) \n',t), the 

variation 5^(t) in Eq. (Bl) is performed separately with respect to the expansion coefficients 
{Cfi(t)} and with respect to the orbitals {(f) k (r,t)}. 

We begin as in the main text with the coefficients. From the structure of \^(t)) as a sum 
of permanents one obviously has 

<*(*)! 



(n;t\ 



(B2) 



dCt(t) 

Hence, making use of the orthonormality relation of the permanents (ft; t \ n'\ t) = S^h 1 one 
straightforwardly finds: 
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H(t)C(t) 
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(B3) 



which is precisely the result depicted in Eq. (10). 

To proceed for the orbitals, we need a few identities and relations. Starting from the basic 
definitions *(r) = ^2 k b k (t)4> k (r,t) and b k {t) = J 4>* k (r, t)^(r)dr connecting the annihilation 
and field operators one readily has, 

Sb k (t) 



mr,t) 



= S kq ^(r). 



(B4) 



With the help of (B4) we can take the variation of a single permanent (n;t\, 

5 (n; t\ 



and from it of the 'bra' ($(t)\ itself, 

5(*(t)\ 



(n;t\bl(t)*(r), 



(B5) 



(B6) 



On the 'kef side of Eq. (Bl) we keep in mind that, 
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dt 
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(B7) 
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Utilizing Eqs. (B3), (B6) and (B7) and the usual commutation relations between the bosonic 
annihilation and creation operators, we evaluate separately the contributions coming from 
the one-body h — and two-body W operators. Combining the results one finds: 



M) 



) = o 



M [ ( d\ M 1 M I 
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l di> bs 



* ) \4>s) -(B8) 



To show that (B8) is exactly Eq. (22) we evaluate the commutator on the right-hand-side 
of the former, 



d ' 

H ~ i W bi 
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(B9) 



31 j 



from which we readily get that the right- hand- side of (B8) is nothing but (sum over) Hkj(t) 
of Eq. (22), i.e., that 



b{ 



d 



(BIO) 



From this point on the two derivations coincide, leading as expected to identical equations- 
of-motion of MCTDHB. 

We conclude the appendix by pointing out an additional connection between the two 
variational formulations used in this work to derive the MCTDHB equations-of-motion. 
Eqs. (B8), (B9) and (BIO) hint towards another role played by the Lagrange multipliers 
Hkj{t) employed in the main text. Since in the Dirac-Frenkel formulation the variation is 
taken before matrix elements are computed, whereas in the Lagrangian formulation the 
situation reverses, there would have been 'missing' terms in the equations-of-motion derived 
by the Lagrangian formulation. The introduction of the Lagrange multipliers in Eq. (6) 
exactly 'compensates' for these terms. 
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X 



FIG. 1: (Color online) Many-body time evolution of two condensates made of N = 100 (in green) 
and N = 1000 (in blue) bosons in a double-well potential with corresponding interaction strengths 
= jfzz computed by the MCTDHB(2) theory. For comparison, the corresponding solution 
computed by the time-dependent Gross-Pitaevskii equation is also displayed (in red). Shown are 
three time snapshots of the density p(x,t), normalized to 1. At t = all densities are identical on 
the scale show. To guide the eye, the double-well trap potential in which the condensates evolve is 
also illustrated (in black). The quantities shown are dimensionless. See text for more details. 
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FIG. 2: (Color online) Natural occupation numbers of two condensates made of N = 100 (in green) 
and N = 1000 (in blue) bosons in a double- well potential (same parameters as in Fig. 1). Shown 
in percents is the largest occupation p\ (p2 = 100% — p\\ P2 is not shown). As time passes, the 
natural occupation numbers change substantially. For comparison, the corresponding fixed 100% 
occupation of the time-dependent Gross-Pitaevskii theory is also displayed (in red). The quantities 
shown are dimensionless. See text for more details. 
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